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Monte Carlo simulation of the nonadditive restricted primitive model of ionic fluids: 

Phase diagram and clustering 

Riccardo Fantono 

Dipartimento di Scienze dei Materiali e Nanosistemi, Universita Ca' Foscari Venezia, 

Calle Larga S. MaHa DD2137, 1-30123 Venezia, Italy 

Giorgio PastorCl 

Dipartimento di Fisica dell' Universita and lOM-CNR, Strada Costiera 11, 34151 Trieste, Italy 

(Dated: April 26, 2013) 

We report an accurate Monte Carlo calculation of the phase diagram and clustering properties 
of the restricted primitive model with non-additive hard-sphere diameters. At high density the 
positively non-additive fluid shows more clustering than the additive model and the negatively 
non-additive fluid shows less clustering than the additive model, at low density the reverse scenario 
appears. A negative nonadditivity tends to favor the formation of neutrally charged clusters starting 
from the dipole. A positive nonadditivity favors the pairing of like ions at high density. The critical 
point of the gas-liquid phase transition moves at higher temperatures and higher densities for a 
negative nonadditivity and at lower temperatures and lower densities for a positive nonadditivity. 
The law of corresponding states does not seem to hold strictly. Our results can be used to interpret 
recent experimental works on room temperature ionic liquids. 
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I. INTRODUCTION 

Ionic soft matter [l| is a class of conventional con- 
densed soft matter [^ |3| with prevailing contribution 
from electrostatics, in many cases crucially affecting its 
physical properties. Among the most popular representa- 
tives of such a class of materials are natural and synthetic 
saline environments, like aqueous and non-aqueous elec- 
trolyte solutions and molten salts as well as a variety of 
polyelectrolytes and colloidal suspensions. Equally well 
known are many biologically important proteins. 

The simplest theoretical model for ionic colloidal sus- 
pensions is the Restricted Primitive Model (RPM) Q, a 
binary mixture of uniformly charged Hard-Spheres (HS) 
of diameter a: two species of opposite charge ±q and 
equal concentrations to ensure charge neutrality, moving 
in a medium of fixed dielectric constant e. The phase 
diagram properties of this model have been widely stud- 
ied both through analytical theories |5l-[l7| and within 
computer experiments starting from the seminal works 
of Larsen [l^l and Vel' Yaminov jig . [20j , followed by the 
pioneering Gibbs Ensemble Monte Carlo (GEMC) cal- 
culation of Panagiotopoulos |2l| and by other numerical 
simulations [22, ,2£] . The more general primitive model 
(PM) with asymmetry in ion charge (30[, in ion size |3ll - 
|33| and in both [S^, ISJ] has also been studied. 

From these studies emerged how, in the vapor phase, 
an important role is played by association and cluster- 
ing. In an old paper [35| one of us studied a modified 
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RPM fluid where one allows for size nonadditivity parti- 
cle diameters. Controlling the nonadditivity, it was sug- 
gested through the use of integral equation theories, that 
such a fluid might have a complex behavior due to the 
possible competition between clustering tendence due to 
the Coulomb interaction and dcmixing tendence due to 
entropic advantage driven by the nonadditivity. Thus, 
the nonadditivity of the hard-sphere diameters does not 
destroy the simplifying symmetry of the model but it 
enriches the properties of the pure RPM model making 
it a paradigm for the self-assembly of isotropic particles 
and a challenge to present day theories of fluids. In real 
systems, the degree of nonadditivity might be directly 
related to the anion-cation contact-pairing affinity |36j 
which in turn may be mediated by the solvent. 

It is the purpose of this paper to reconsider such a 
model fluid from the point of view of accurate numer- 
ical experiments. In particular, we want to study the 
clustering properties of the fluid outside of the gas-liquid 
coexistence region. To this aim we first determined the 
gas-liquid coexistence curve through the Gibbs ensem- 
ble method after having studied semi-quantitatively how 
the coexistence region changes with the nonadditivity 
through a density distribution analysis in the canonical 
ensemble. This way we could be sure that our cluster 
analysis falls outside the coexistence region in all the 
cases studied. Clustering turns out to be greatly affected 
by the nonadditivity parameter. The most striking effect 
being the prevalence of neutrally charged clusters made 
up of an even number of particles in the negatively non- 
additive fiuid. When the nonadditivity allows complete 
overlap of the two species of particles the formation of 
a fluid of neutral hard-spheres of half the density is ex- 



pected and our simulation results clearly show this be- 
havior. On the other hand, for a positive nonadditivity, 
it is known that the neutral HS mixture tends to demix 
the two species and the demixing critical density low- 
ers as the nonadditivity increases |37[. We expect this 



property of the neutral system to have some interesting 
effect on the clustering properties of the charged fluid 
since demixing cannot occur in a binary charged system: 
The frustrated tendency to segregation of like particles 
and the reduced space available to the ions favors pairing 
of like ions and percolating clusters at high densities. To 
the best of our knowledge this is the first time that such 
a model fluid is studied with numerical simulations. Pre- 
liminary results from our analysis have been presented 
in a letter |38| , here we extend that analysis and present 
for the first time the gas-liquid binodal of the fluid as a 
function of the nonadditivity parameter. 

We think that the model fluid considered in this pa- 
per may be realized experimentally through a colloid-star 
polymer mixture where both species are charged |39[ or 
by room temperature ionic liquids |40l443J as discussed in 
Section lilIC2l In particular in the latter systems liquid- 
liquid binodals shifted above and below the one of the 
pure RPM are observed depending on the kind of sol- 
vent used. If on the one hand this can be ascribed to 
the different dielectric constant of the solvent 141[, on 
the other hand it is clear that, depending on the kind 
of solvent, the anion-cation contact-pairing affinity may 
vary [Sa] and thus the different experimental ionic liq- 
uids should be more correctly described by comparison 
not just with the pure RPM but with the more realistic 
primitive model with the addition of either a positive or 
negative size nonadditivity. 

The paper is organized as follows: in Section |ll] the 
model for the fluid we want to study is described, in Sec- 
tion lllll thc results from the numerical experiments arc re- 
ported. These are divided in a cluster analysis in Section 
nil Al in an analysis of the radial distribution function 
and structure factor in Section IIIIBl and in an analysis 
of the gas-liquid coexistence in Section IIII CI Theoreti- 
cal remarks on the clustering properties are presented in 
Section ITVl and Section |V] is for final remarks. 



II. THE MODEL 



The model fluid we want to study is the restricted 
primitive model (RPM) of non-additive hard-spheres 
(NAHS). The RPM consists of N/2 uniformly charged 
hard-spheres of species 1 of diameter a carrying a total 
charge +q each and N/2 uniformly charged hard-spheres 
of species 2 of the same diameter carrying a total charge 
—q each. The spheres are moving in a dielectric contin- 
uum of dielectric constant e independent on the thermo- 
dynamic state. The interaction between an ion of species 



i and one of species j a distance r apart is given by 
f +00 r < Oi 



/30^j(r) 



'%2 
7ij 



g'gj r>a ' *'J' = 1'2, (2.1) 
Kb! er 



where /3 — I/fc^T with T the absolute temperature and 
ks the Boltzmann's constant, qi the charge of an ion of 
species i. The ions form a mixture of non- additive hard- 
spheres, i.e. 



a 1 = ] 

(t(I + A) i^j 
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1,2 



(2.2) 



with the nonadditivity parameter A > — 1 . A thermody- 
namic state is completely specified by the reduced density 
p* = pa^ = N(j^/V, where V is the volume contain- 
ing the fluid, the reduced temperature T* = ksTea/q'^ 
(g^/(ecr) is taken as unit of energy), and the nonaddi- 
tivity parameter A. We will call xi — p^/p = 1/2 and 
X2 = P+/p = 1/2 the anions and cations molar concen- 
trations respectively. 



III. RESULTS 

In Fig. [T] we show the phase diagram of the pure RPM 
fluid, A = 0, as obtained from the Gibbs ensemble Monte 
Carlo method by Orkoulas et al. [2l] and by us (see 
Section FlIIC 2p . The thermodynamic points where we 
probe the fluid with our NVT Monte Carlo simulations 
are also shown as filled squares. 
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FIG. 1. (Color online) Diagram showing the gas-liquid coex- 
istence curve of the RPM fluid from the Gibbs ensemble MC 
data of Ref. [2J] (filled circles) and ours (open circles), the 
triangle being the critical point, and the points (filled squares) 
of the phase diagram where we run our NVT MC simulations. 

In our canonical NVT Monte Carlo (MC) we study 
the fluid in a simulation box of volume V ~ L^ with pe- 
riodic boundary conditions. The long range of the 1/r 
interaction is accounted for using an Ewald sum for the 



interacting energy in the periodic system [441 . The inter- 
action energy per unit box for e = 1 is calculated as 



,-^ .^-^ erfc(K|r^,, + Ln 



/^<z^ n 



Tf^^ + Ln\ 



2n er^kl2.f 



^(E^v) - (3.1) 
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where a Roman index with a Greek sub-index denotes 
the species of the particle labeled by the Greek sub- 
index, r^i, = r^ — r^ with r^ the position of particle 
M) Pk = X]u 9^^ *'^'""i "^rfc denotes the complemen- 
tary error function, n = (77,3;, riy, n^) with nx,ny,nk = 
0, ±1, ±2, . . ., and k = (2-e j V){nx, Uy, Uz) are reciprocal 
lattice vectors. The parameter k which governs the rate 
of convergence of the real space and reciprocal space con- 
tribution to the energy is taken to be k ~ 5/L. With this 
value of K, the real space contribution can be restricted 
to the first term n = only. The reciprocal space term 
includes all k vectors such that ni + nf, + n1 < 27. The 
last term in Eq. (|3.ip is zero for the RPM but it is impor- 
tant in the Gibbs ensemble simulation where a particle 
exchange between the two boxes can produce systems 
where there is an unequal number of positive and nega- 
tive charges. Our choice for the interaction energy takes 
into account the fact that each charge has a uniform back- 
ground of neutralizing opposite charge density. 

In our NVT MC simulations we used N = 100 (ex- 
cept for the test of the size dependence of the clustering 
analysis where we considered up to 5000 particles), the 
acceptance ratio is kept, on average, close to 50% after 
a preliminary adjustment of the maximum particle dis- 
placement. We start from a simple cubic configuration 
of two crystals one made of species 1 and one made of 
species 2 juxtaposed in order to avoid overlaps at high 
densities. We need around 10^ MC steps (MCS) in order 
to equilibrate the samples and 10^ MCS/particle for the 
statistics. 



A. Cluster analysis 

During the simulation we perform a cluster analysis. 
After each 100 MCS we determine the number iV„ of 
clusters made of n particles, so that X^n '^^" ~ ^- ^'^ 



45| that a group of ions forms a cluster if the dis- 



tance r, calculated using periodic boundary conditions, 
between a particle of species i of the group and at least 
one other particle of species j is less then some fixed 
value, i.e. r < (Tij + S'^cr where 5'^ is a paramete r l46i . 
In all our simulations we choose J'^ = 0.1 (in Ref. [23 a 
detailed study of the sensitivity of the clustering proper- 
tics on this parameter is carried out for the pure RPM 
fiuid). Then we take the average of these numbers (iV„). 



Note that Q„ = n{Nn)/N gives the probability that a 
particle belongs to a cluster of size n. To establish a cri- 
terion for percolation we first find the clusters without 
employing periodic boundary conditions to calculate the 
distances, then we check whether, amongst the particles 
of any of these clusters, there are two which satisfy the 
cluster condition calculating the distances using periodic 
boundary conditions. Whenever we find one such cluster 
the cluster is percolating. 

In Fig. [5] we show results of such analysis for the fluid 
at a temperature T* = 0.1 well above the critical tem- 
perature, T* w 0.05, of the pure RPM [H, [H, HI] . In 
the insets wc show a magnification of the region around 
n = 1 from which the degree of dissociation [ij, ll3| 
a ~ {Ni)/N can be read-off. In the figure we plot the 
cluster concentrations {Nn)/N as a function of the num- 
ber of particles n in the cluster. We plot n from n = 1 
(isolated ions) up to tt, = A'' (in this case all the parti- 
cles of the fluid form one big percolating cluster). At 
p* == 0.45 both the pure RPM and the A = +0.3 fluid 
form percolating clusters. Lowering the density we first 
reach a state, at p* = 0.3, where the negative nonaddi- 
tivity gives the same clustering of RPM and the positive 
nonadditivity gives bigger clustering (still with percolat- 
ing clusters), then a state, at p* =0.1, where the positive 
nonadditivity gives the same clustering of RPM and the 
negative nonadditivity a bigger one, and finally a state 
p* = 0.01, 0.001, at low densities where a negative nonad- 
ditivity increases the clustering over the RPM fiuid and 
a positive nonadditivity diminishes it. Generally, at high 
densities we find percolating clusters in the fluids whereas 
these disappear at low densities even at a positive non- 
additivity. Summarizing, in agreement with Ref. [35|, 
we find, for the fixed values of |A|, that: at high density 
and positive A we have more clustering than in the ad- 
ditive model since there is a smaller effective volume for 
particle, at high density and negative A we have less clus- 
tering than in the additive model because there is more 
effective volume for the particles, at low density and pos- 
itive A we have less clustering than in the additive model 
due to the competition between tendency to demixing in 
the corresponding neutral mixture and tendency to local 
electroneutrality of the Coulombic systems, at low den- 
sities and negative A we have more clustering than in 
the additive model because neutral clusters are favored 
as shown in the next Section. We conclude that at high 
temperature and high density the negative nonadditiv- 
ity gives lower clustering than RPM, lowering the tem- 
perature at constant density or lowering the density at 
constant temperature it gradually tends to gives higher 
clustering than RPM. On the contrary, at low density the 
positive nonadditivity gives lower clustering than RPM, 
increasing the density it gradually tends to give larger 
clustering than RPM. 

We determine the size dependence of the curves shown 
in Fig. [21 and see that when we have no percolating clus- 
ters, for example the data at T* = 0.1, p* = 0.3, A = 
0, —0.3, the curves were unaffected by a choice of a higher 



number of particles, while when we have percolating clus- 
ters, for example the data at T* = 0.1, p* = 0.3, A = 0.3, 
the curve (n, {Nn)/N) changes with N. In these lat- 
ter cases we find that a common curve is given by 
(x, {Nx)/N) with X = n/N e [0, 1]. Then, in order to sat- 
isfy the normalization condition, 1 = X^n "^((^n)/^) ~ 
J dxxN'^{{Nx)/N), we must have for two different sizes 
A^' andiV" that {{N^)/N')/{{N^)/N") w {N"/N'f. We 
have no general recipe to when the former behavior is to 
be expected over the latter. We can only say that the first 
behavior is generally observed when we do not have per- 
colating clusters whereas the second is present when we 
have percolating clusters. In Section ITVl we show that the 
size independent curves that we find when there are no 
percolating clusters can be fitted by {Nn)/N — a"-n''"-/nl 
(see Eq. (|4.4p with 2™*''" obtained from an ideal cluster 
approximation) with a and b a positive fitting parame- 
ter. In Table |T] we show the fitting parameters a and b 
corresponding to the simulated cases. 



TABLE I. Fitting parameters a, b in the least square fit 
{Nn)/N = a"n'"'/n\ for the simulation results of Fig. [2]with- 
out percolating clusters (and with tire exclusion of the non 
smooth data at p* = 0.001). The reduced x^ was around 
0.5 with greater error approaching n = 1. Also shown is the 
number of particles Umax in the biggest cluster formed in each 
simulation. 
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^max 


0.45 


-0.3 


0.220(3) 


1.074(4) 


64 


0.3 





0.197(4) 


1.084(6) 


45 


0.3 


-0.3 


0.204(3) 


1.069(5) 


43 


0.2 





0.206(7) 


1.00(1) 


23 


0.2 


-1-0.3 


0.200(4) 


1.083(5) 


45 


0.2 


-0.3 


0.204(7) 


1.04(1) 


31 


0.1 





0.22(2) 


0.86(3) 


15 


0.1 


-1-0.3 


0.16(1) 


1.01(4) 


19 


0.1 


-0.3 


0.15(1) 


1.11(2) 


29 


0.01 





0.41(7) 


0.1(1) 


8 


0.01 


-1-0.3 


0.36(8) 


0.0(2) 


7 


0.01 


-0.3 


0.23(4) 


0.72(7) 


12 



In Fig. |3]we show the clustering analysis at the ther- 
modynamic state below the critical temperature of RPM 
T* = 0.04 in the gas phase, p* = 5 x 10^^, and in the 
liquid phase, p* = 0.45. We see how in the gas phase 
only the first few clusters are present in agreement with 
similar results found in Ref. [25| and for a negative non- 
additivity the dipoles are clearly the preferred kind of 
clusters with the smallest degree of dissociation amongst 
the three fluids considered. In the liquid phase all three 
fluids have percolating clusters. 

In Fig. |4]we show the clustering analysis for the fluid 
with A approaching -1 at T* = 0.1 and p* = 0.45. We 
see how letting A approach —1 this stabilizes the neu- 
trally charged clusters and lowers the degree of dissocia- 
tion. The first stable cluster is the dipole: the "overlap" 
of a positive and a negative sphere. These are dipoles of 



moment qri2 with ri2 < cr(l -|- A -|- 5'^) which may lack 
a gas- liquid criticality [43]. We clearly have a transition 
from a conducting to an insulating phase as A goes from 
to -1. 

In Fig. [5] we show a snapshot of the equilibrated fluid 
at T* = 0.1, p* = 0.45 and A = -0.9 from which one 
can sec the formation of the dipoles. We expect that in 
the limiting case of A — —1 the fluid we obtain is well 
reproduced by hard-spheres at half the density. This is 
confirmed by a comparison of the like radial distribu- 
tion functions with the one of the hard-spheres even if 
the A = — 1 fluid simulation rapidly slows down into 
the frozen configuration of the overlapping anions and 
cations. In order to overcome this problem one should 
alternate single particle moves to cluster moves where 
one moves the center of mass of the neutrally charged 
pairs. 



B. Radial distribution function and structure 
factor 

In Figs. [S][5] we show the partial radial distribution 
functions (RDF), g,,{r) = {Y.'^,J{^+^i~^l))N /{px,x,), 
where r* denote the position of particle p of species i 
and the prime over the sum indicates that the terms p = 
V when i = j are omitted, and the total RDF, gtot = 
Si j=i9ij^i^j^ of the three fluids A ~ 0,±0.3 at the 
thermodynamic states T* = 0.1, p* — 0.01, 0.1 and T* ~ 
0.04, p* = 0.45. Of course, the restrictions xi — X2 and 
fii = 1722 imply 5ii = (?22- In the simulations we use 
N = 100. 

From Fig. [S] we see how the contact value of the like 
RDF in the A = —0.3 case is higher than in the additive 
case and in the A = -1-0.3 case is lower than in the addi- 
tive case. The contact value of the unlike RDF is high- 
est for negative nonadditivity indicating the tendency to 
form cation-anion pairs. 

From Fig. [7] we see again the same behaviors of the 
contact values of the like and unlike RDF. In the neg- 
atively additive case we begin to see an alternation of 
distribution of oppositely charged shells of ions around a 
reference ion. 

From Fig. [5] we see how at this high density the con- 
tact value of the like RDF is highest in the A = —0.3 
case but in the A = -1-0.3 case is still higher than in the 
additive case. At A = —0.3 we see clearly the formation 
of a second peak in the unlike RDF around 2 + A and 
the expected alternation between the peaks of the like 
RDF with the ones of the unlike RDF also present in the 
additive case. This alternation is not present in the pos- 
itively non-additive case indicating now the tendency of 
like particles to cluster on a microscopic scale: like parti- 
cles penetrate inside the shell of unlike particles around 
a given reference ion. The contact value of the unlike 
RDF is highest for negative nonadditivity indicating the 
tendency to form cation-anion pairs. 

In Fig. [5] we show the Bhatia-Thornton [4811 structure 
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FIG. 2. (Color online) Clustering properties of the fluid at T* — O.f at various values of nonadditivity. Nn are the number of 
clusters made of n particles. In the MC simulations we used A'' = 100 particles and a number of MCS= 1 x 10 . The panels 
are ordered (left to right, top to bottom) in order of decreasing density p* — 0.45,0.3,0.2,0.1,0.01, and 0.001 for panels (a), 
(b), (c), (d), (e), and (f) respectively. The insets allows to read-off the degree of dissociation. 



factors Snn = [Sii + S22 + 2S'i2]/2 and Sqq = [6*11 + 
S22 — 25i2]/2 where Sij{k) = {pi_p'_]J / {N ^/xixj) are the 
partial structure factors and pj^. — ^ exp(— ik-rM is the 
Fourier transform of the microscopic density of particles 
of species i. In the figure we chose the same thermo- 
dynamic state and nonadditivity considered in Fig. 1 of 
Ref. [Sa] . The positive nonadditivity case has percolating 
clusters. From the figure we see that the charge-charge 
structure factor Sqq tends to zero at fc = a conse- 



quence of electroneutrality in charged systems [J] which 
suppresses long-wavelength fluctuations. In order to en- 
force this condition the structure factor needs to develop 
a peak at small k which reflects an essentially alternating 
distribution of oppositely charged shell of ions around a 
reference ion. This type of short range order is an indi- 
cation of the tendency to cluster. From the figure we see 
that at high density the positive non-additive fluid tends 
to cluster more than the additive fluid and the negative 
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FIG. 3. (Color online) Clustering properties of the fluid at 
T* = 0.04 and p* = 0.45,5 x fO~^ for panels (a) and (b) 
respectively, at various values of nonadditivity. Nn are the 
number of clusters made of n particles. In the MC simulations 
we used A^ = 100 particles and a number of MCS= 1 x 10^. 
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FIG. 4. (Color online) Clustering properties of the fluid at 
T* = 0.1 and p* — 0.45 at various values of negative non- 
additivity approaching —1. Nn are the number of clusters 
made of n particles. In the MC simulations we used A'' = 100 
particles and a number of MCS= 5 x 10^. 




FIG. 5. (Color onhne) Snapshot of the fluid at T* = 0.1, p* = 
0.45 and A = —0.9 showing the formation of the dipoles. 



non-additive fluid tends to eluster less than the addi- 
tive fluid, in agreement with the results presented in the 
previous Section. With regard to the number-number 
structure factor Snn we see that as the nonadditivity 
decreases the isothermal compressibility Snn{0) (see the 
appendix of Ref. [4^ |4^ ) increases and the short range 
order is reduced. 

In Table |ll] we report the excess internal energy per 
particle W^^/N ~ ea{U)/Nq^, the compressibility fac- 
tor Z = pP/p, and the total clusters concentration 
J2ni-^n)/N for the cases simulated. The compressibil- 
ity factor is calculated according to the virial theorem 



U'' 



3NT* 



+ ^[.9ii(^) + (l + A)3.gi2Kl + A))] . (3.2) 



If the clusters do not interact, as in the independent clus- 
ter model (ICM) of GiUan 0, one should have Z^'^^ = 
Tlini-^n) /N . From Table ITTl we can see how this condition 
is never satisfied in the cases considered. 



C. Gas-liquid coexistence 

An important question we try to answer is how the gas- 
liquid coexistence curve of the pure RPM fluid changes 
upon switching on of the nonadditivity parameter. To 
this aim we first perform a density distribution analy- 
sis within the NVT ensemble which allows us to easily 
extract a semi-quantitative result and then we use the 
Gibbs ensemble technique for a careful quantitative de- 
termination of the binodals. 
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FIG. 6. (Color online) Partial and total RDF in the simulations at p* = 0.01, T* = 0.1, and A = (panel (a)), A = 0.3 
(panel (b)), A = —0.3 (panel (c)). The reduced excess internal energy per particle of the fluid was in the three cases 
f/™/iV = -0.3924(1), -0.29120(7), -0.6339(1) respectively. 



1. The density distribution approach 

Sufficiently close to the critical point we determine how 
semi-quantitatively the behavior of the gas-liquid coexis- 
tence region changes by switching on a negative or a pos- 
itive nonadditivity. To this aim we divide the simulation 
box into 771"^ cubes of side Lc = Ljm and register, as the 
run progresses, the density inside each cell pi = Mi/ L^, 
where A/i is the number of particles inside the ith cell so 

3 

that YlT=i-^i ~ ^- Then we calculate the density dis- 
tribution function [sfll Pm{p) = Z^ili Pm{pi)/m^, where 
Pm{pi) is the distribution function for the ith cell, with 
/ Pm {p) dp = 1. Above the critical temperature the den- 
sity probability distribution function can be described by 
a Gaussian centered at the simulation density whereas 
below it becomes bimodal with two peaks one centered 
at the gas density and one at the liquid density. 

We start from an initial configuration of particles of 
random species placed on a simple cubic lattice. We 
equilibrate (melt) the fluid for 10*^ MCS/particle. We 
then sample the distribution function every 10 MCS. To 
allow the particles to diffuse out of the cells we choose the 
subdivision of the simulation box in cells with a random 



displacement r = {rx,ry,rz) with r,j;,ry,rz € [0,i]. This 
procedure turned out to greatly enhance the efficiency of 
the determination of the cell density distribution. And 
we measure the distribution function on runs of 1 x 10^ 
MCS/particle. 

Choosing m = 2 and A^ = 100 the results for the fluid 
at a temperature T* = 0.025 above the triple point of 
the RPM [5l|, a density p* = 0.2 well within the coex- 
istence region of the pure RPM fluid, and A = 0, ±1? 
with V = 10-\l0-2,5 X 10-2 are shown in Fig. [TUl 
In this case the minimum density that can be registered 
is 1/Ll = 0.2 X 8/100 = 0.016. We see that the pure 
RPM fluid shows a density distribution function with two 
peaks: The first one, which lies below the minimum den- 
sity (and is not visible in our data) , at approximately the 
low density of the gas phase and the second one at ap- 
proximately the high density of the liquid phase around 
a reduced density of 0.3. At 2? = 10"^ the positions of 
the peaks are roughly the same to the pure RPM. At 
T) = 5 X 10^2 the density of the hquid peak in the neg- 
atively non-additive fiuid is higher than the one of the 
pure RPM whereas the positively non-additive fluid has 
a gas peak, now visible, at higher density than the pure 
RPM and a liquid peak at lower density than the pure 
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FIG. 7. (Color online) Same as Fig. [6]at p* = 0.1 and T* = 0.1. The reduced excess internal energy per particle of the fluid 
was in the three cases U""^ jN = -0.50589(8), -0.41208(6), -0.7179(1) respectively. 



RPM. At P = 10^^ this separation tends to increase: In 
the positively non-additive model the eritieal tempera- 
ture is too close to 0.025 and the bimodal is degenerate 
into a curve with a single peak centered on the simula- 
tion density 0.2 whereas in the negatively non-additive 
fluid the liquid peak is changed into a broad tail extend- 
ing up to a density of 0.8. This findings suggest that at 
a given temperature the width of the coexistence region, 
relative to the one of the pure RPM, tends to increase 
for the negatively non-additive model and to decrease for 
the positively non-additive model. This result is made 
more clear and precise in the following Section where we 
present our Gibbs Ensemble Monte Carlo calculation. 



2. Gibbs ensemble analysis 

In order to quantitatively determine the gas-liquid co- 
existence line of our fluid we use the Gibbs ensemble MC 
(GEMC) technique [s^], starting from the pure RPM one 
and gradually switching on the nonadditivity. Here, we 
are not interested in the behavior really close to the crit- 
ical point but rather in the shape of the binodal curve 
and how it moves as a function of A. 

The GEMC method of Panagiotopoulos is now widely 



adopted as a standard method for calculating phase equi- 
libria from molecular simulations. According to this 
method, the simulation is performed in two boxes con- 
taining the coexisting phases. Equilibration in each 
phase is guaranteed by moving particles. Equality of 
pressures is satisfied in a statistical sense by expanding 
the volume of one of the boxes and contracting the vol- 
ume of the other. Chemical potentials are equalized by 
transferring particles from one box to the other. Like 
the first simulations for the RPM performed by Pana- 
giotopoulos [21j we use single ion transfer by introducing 
a background charge density to assure charge neutrality 
at all times during the run. This way, the system re- 
mains overall neutral, but the modified model is similar 
to a two component plasma and in a strict sense differ- 
ent from the original RPM model which assumes a zero 
charge density for the background. To overcome the elec- 
troneutrality problem Orkoulas et al. |2j| considered pair 
transfers. 

In the GEMC run we have at each step a probabil- 
ity ap/{ap + a^ + Us) for a particle random displacement, 
CLv/iap+ay+as) for a volume change, and as/{ap+ay+as) 
for a particle swap move between the gas and the liq- 
uid box. We generally choose Up = f,a„ = 1/10, and 
as = 1. The maximum particle displacement is kept 
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equal to Li/1000 where Li is the side of the ith box with 
i = 1,2. Regarding the volume changes, following Ref. 
[52| we perform a random walk in ln[V^i/V2], with Vi the 
volume of the ith box choosing a maximum volume dis- 
placement of 1% — 10%. Volume moves are computation- 



ally the cheapest since the energy scales with the length 
of the box with inverse proportionality. We generally use 
a total number of A^ = 100 particles except close to the 
critical point where it proves necessary to increase the 
number of particles in order to avoid large fluctuations 
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in the two densities. We use 10 — 40 million MCS for the 
equilibration and 100 — 200 million MCS for the produc- 
tion. 

The results are summarized in Table IIIII and Fig. [TT] 
Note that since we get the same coexistence curve of Ork- 
oulas et al. 2% for the pure RPM, as Fig. [T] clearly 
shows, we consider as equivalent, at the present level 
of accuracy, our procedure, employing single neutral- 
ized particle transfers, and the one of Orkoulas, where 
pair particle transfers between the two boxes are used. 
This can be justified by observing that the fluctuations 
of charge in the various statistical physics ensembles are 
expected to decay to zero with the system size and we 
empirically find that in our case they are already practi- 
cally irrelevant. The only relevant difference we observe 
with respect to the calculation of Orkoulas is the fact 
that in our case there is a much more considerable emp- 
tying of the gas box at low temperatures which may have 
some effect in the point at the lowest temperature. We 
do not carry out a systematic study of possible system 
size dependence of the results but, for the pure RPM, we 
repeat the calculation at T* = 0.045 and 0.0475 for two 
different system sizes with the largest being iV = 370. 
The comparison suggests that the critical point tends to 
shift slightly at higher temperatures upon a system size 



increase but far away from the critical point the coex- 
istence curve is not affected appreciably by the system 
size. However, we stress that an accurate study of crit- 
ical properties of the present model is beyond the scope 
of this work. 

From the figure we can see clearly the trend: a pos- 
itive nonadditivity tends to lower the critical tempera- 
ture whereas a negative one tends to push the binodal 
to higher temperatures. This is in agreement with the 
findings from the density distribution analysis previously 
presented. It is well known that RPM condensation is 
almost identical to that of charged hard dumbbells, un- 
derlining the fact that the vapor is essentially already 
fully associated into dimers and higher neutral clusters, 
and that the liquid structure and thermodynamics are 
only weakly perturbed by fusing ions together. Hence, 
if one imagines cooling down on the critical isochore, we 
can say that the critical point is reached when ion associ- 
ation is complete and then it becomes convenient for the 
system to phase separate. With positive nonadditivity, 
ion association is less favorable and the critical tempera- 
ture must go down (association is complete only at lower 
temperatures); with negative nonadditivity, ion associa- 
tion is more favorable, and the critical temperature must 
go up (association is complete already at higher temper- 
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TABLE II. Excess internal energy per particle C/'^" = 
ecr(bl)/q^, compressibility factor Z — pP/p, and total clus- 
ters concentration '^„{N„)/N for the cases simulated. 
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In order to determine the critical point (T^*,p*) we 
empirically fit the binodals using the "diameters" (p* + 
p1)/2 equation [53J 



P*g + P*l 



^|7^* - T;| + C\T* - T*\^^' 



+D\T* -T*\^-"" , (3.3) 

and the form of the Wegner expansion [SJ, [SJj for the 
width of the coexistence curve 



P;* - p*g = B\T* - T*\f^' + Bi\T* - T2f'+^' 



(3.4) 



where A,C,D, and B,Bi,B2 are coefficients which we 
take as fitting parameters as well as pl,T*. We stress 
that our data do not extend sufficiently close to the crit- 
ical region to allow quantitative estimates of critical ex- 
ponents and non universal quantities, still we used the 
above functional forms as convenient fitting formulae, 
able to capture the typical flatness of the fluid coexis- 
tence curves ^. The pure RPM is believed [H [H-IE^I 
to belong to the three-dimensional Ising universality class 
so we choose /3j ~ 0.325, a/ = 0.11, and A/ = 0.51. 
We are then able to flt the pure RPM case, A = 0, 
for which wc find the critical point at p* = 0.0319 
and T* = 0.0476, the RPM with positive nonadditiv- 
ity, A = -1-0.1, for which the critical point is found at 
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FIG. 11. (Color online) Binodals obtained with the GEMC 
simulations. The dashed lines are the result of the extrapo- 
lation described in the text. The filled circles are the critical 
points. 



p* = 0.0275, T* = 0.0432, and the RPM with nega- 
tive nonadditivity, A = —0.1, for which p* = 0.0495, 
T* ~ 0.0526. We stress that these numbers, in partic- 
ular the values of critical densities should be considered 
more as indicative of the dependence of the critical point 
location on diamater non additivity than as accurate es- 
timates. 

We believe that our results can be relevant for the 
interpretation of experimental work on the phase dia- 
grams of room temperature ionic liquids [42| like the 
phosphonium halogenide in alkanes solvents and 1-hexyl 
3-methyl imidazolium tetrafluoro borate (C6niimBF4) in 
alcohols and water. The degree of nonadditivity seems di- 
rectly related to the anion-cation contact-pairing affinity 
[3a |. The salts in the (hydrocarbon) solution dissociate 
in cations (the phosphonium) and anions (the halogen 
atoms). The contact affinity between anions and cations 
is mediated by the solvent and different solvents produce 
different affinities. As a consequence, in the experimental 
work of Ref. [i^l they observe liquid-liquid coexistence 
curves which, depending on the kind of solvent used in 
the ionic liquid mixture, can be above (C6mimBF4 in 
alcohols and water) the one of the pure RPM theoreti- 
cal model or below (phosphonium halogenide in alkanes) 
in reduced units. Moreover, when plotted into a cor- 
responding state representation all the experimental bin- 
odals seem to collapse on a same curve even if this occurs 
very close to the critical point. We then try to see if the 
law of corresponding states holds or not for our fluid and 
we find that far from the critical point it is not strictly 
satisfied, as shown by Fig. [T^l Interestingly enough, 
a plot of the RDF between corresponding states shows 
an almost complete overlap of the three curves upon a 



12 



TABLE III. Phase coexistence properties for the pure RPM (A — 0) and the non-additive RPM (A 7^ 0). T* is the reduced 
temperature, A^ is the total number of particles in the system for a certain run, Ng is the average number of particles in the 
gas box during the run, fj,* = fiiea/q^ — T* In A^ is the reduced chemical potential of the liquid box (A being the de Broglie 
thermal wavelength), Ui'^ is the total excess internal energy, and p* is the reduced density, of the gas phase i — g and the liquid 
phase i — I. 



A 


T* 


N 


NJN 


~f^i 


-c/r/^ 


~ur/N 


Pi 


Pi 













0.0475 
0.0475 

0.045 

0.045 
0.0425 

0.04 
0.0375 

0.035 


370 
200 
370 
100 
100 
100 
100 
100 


0.51(1) 
0.33(1) 
0.26(1) 
0.27(1) 
0.166(8) 
0.069(5) 
0.036(2) 
0.0020(6) 


0.63(1) 
0.63(1) 
0.69(3) 
0.63(2) 
0.65(1) 
0.73(1) 
0.72(1) 
0.75(2) 


0.547(1) 
0.559(2) 
0.528(4) 
0.537(4) 
0.52(1) 
0.50(2) 
0.4(1) 
0.05(40) 


0.609(1) 
0.604(1) 
0.6400(7) 
0.6393(9) 
0.6576(8) 
0.6745(5) 
0.6835(5) 
0.6938(5) 


5.2(9) X 10"^ 
1.1(3) X 10"^ 
2.3(5) X 10-^ 
3.1(7) X 10"^ 
2.3(4) X 10"-^ 
8(3) X 10"'' 
4(2) X lO-* 
2(20) X 10"^ 


0.11(3) 
0.08(3) 
0.22(5) 
0.22(5) 
0.29(2) 
0.35(3) 
0.38(5) 
0.42(2) 


-0.1 
-0.1 
-0.1 
-0.1 
-0.1 
-0.1 


0.0525 

0.05 

0.0475 

0.045 

0.0425 

0.04 


200 
100 
100 
100 
100 
100 


0.297(9) 
0.37(1) 
0.094(4) 
0.092(5) 
0.031(3) 
0.0034(7) 


0.71(2) 
0.67(1) 
0.71(2) 
0.69(2) 
0.85(3) 
0.83(2) 


0.602(2) 
0.609(4) 
0.562(7) 
0.590(9) 

0.4(1) 
0.08(40) 


0.6844(9) 
0.712(1) 
0.7240(8) 
0.7380(8) 
0.7503(8) 
0.7582(7) 


9(2) X lO"-'' 
4.8(7) X 10"-^ 
1.7(3) X 10"^ 
1.3(2) X 10"^ 
5.6(10) X 10"* 
5(30) X 10"^ 


0.19(5) 
0.33(5) 
0.36(7) 
0.42(5) 
0.46(4) 
0.48(4) 


+0.1 
+0.1 
+0.1 
+0.1 


0.0425 
0.04 

0.0375 
0.035 


100 
100 
100 
100 


0.40(2) 
0.23(1) 
0.068(5) 
0.015(2) 


0.58(1) 
0.62(3) 
0.70(2) 
0.72(2) 


0.493(3) 
0.487(2) 
0.40(4) 
0.17(26) 


0.5620(9) 
0.5877(9) 
0.6068(5) 
0.6157(6) 


5(1) X 10"^ 

1.6(4) X 10"^ 

8(1) X 10"* 

1.7(40) X 10"* 


0.11(3) 
0.19(5) 
0.24(3) 
0.28(2) 



shift by ± A in r, as is show by Fig [131 We think that 
the only visible difference, the contact values of the like 
RDF, is a direct hallmark of the break-up of the corre- 
sponding states, as physical consequence of the existence 
of a third relevant interaction parameter, in addition to 
the unlike hard-sphere diameter and the electric charge. 
While the Cl~ ion and the BF4 anion may reasonably 



charge is identical with the center of mass, the NTF^ an- 
ion is by no means spherical. The NTF^ anion is flexible 
and allows for different conformers. The nitrogen atom 
in the anion is not necessary identical with the center of 
mass and the center of charges [4J]. In these cases in- 
stead of the RPM it is better to choose the PM with ions 
of differing sizes as reference system [3ll - l33| . 
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FIG. 12. (Color online) Corresponding state representation 
of the phase diagram. 



FIG. 13. (Color online) RDF on corresponding states of 
the liquid branch at T' = 0.0375, p* = 0.38, A = 0, T* = 
0.035, p* = 0.28, A = +0.1, and T* = 0.0425, p* = 0.46, A = 
—0.1. In the inset are shown the functions shifted by ±A. 



well be approximated by a sphere so that the center of 
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IV. THEORETICAL REMARKS ON THE 
CLUSTERING 

Under highly diluted conditions [5^ we can approxi- 
mate the fluid as an ideal mixture of cation-anions, an- 
ions, and cations with partial densities /0± = (1 — a)p/2, 
p^ = p^ = Q;p/2, respectively, and for the chemi- 
cal potentials p± = fc_Brin((l - a)pK\K^_/2K), ^_ = 
kBT\\\{apK\/2), and p+ = kBT\n{aph!t/2), where A_ 
and A+ are the de Broglie thermal wavelengths of the 
anions and cations respectively. K is the configurational 
integral of a cation-anion pair 



is: = 47r 



^2gAi 



'"'dr 



cr(l+A) 



(4.1) 



where A^ = cr/T* is the Bjerrum length and Tc is a cutoff 
radius conventionally chosen equal to As/2 correspond- 
ing to the minimum of the integrand. At equilibrium 
p± = /x_|_ + /x_, which implies (1 — a)/oP' ~ Kp/2. Solv- 
ing for a yields 



yi + 2Kp - 1 
Kp 



(4.2) 



An approximate closed form expression for K valid at low 
T* can be obtained by writing for the anion-cation pair 
distance r = (t(1 + l^) + 5r with 5r small. Then a /r « 
1/(1 + A) - ad>/cr2(i + A)2 = 2/(1 + A) - r/a{l + l\f. 
Substituting into Eq. (|4.ip and performing the integral 
with Tc = oo yields 

7^ « 4^^3(1 + A)V/^*(i+^)T* 

x{l + 2(1 + A)T*[1 + (1 + A)r*]} . (4.3) 

In our simulations we are never in this very diluted 
condition and as a consequence we observe the forma- 
tion of clusters of an higher number of particles than 
just the dimers. So to estimate the cluster concentra- 
tions x^ = {Nn}/N, we need a different analysis closer 
in spirit to the one of Tani and Henderson [4^, [S^] . Sim- 
plifying that analysis we can consider as the inter-cluster 
configurational partition function the one of an ideal gas 
of clusters, in reduced units, Zinter ~ (1^/c^)^', where 
Nt = X]n=i ^n i^ *^^ total number of clusters and we 
assume to have only clusters made of up to n'^ particles. 
Then the equations for the equilibrium cluster concen- 
trations x'i are 



A"z!: 



Vp* 



1,2,. 



i = E^ 



(4.4) 
(4.5) 



ri=l 



where z™*''" are the configurational intra-cluster par- 
tition functions in reduced units with 2™*^° = 2 and 
A(= ap* /2) is a Lagrange multiplier. Moreover, ne- 
glecting the excess internal energy of the clusters we can 



approximate z™*""" 

-.3\n — Ion 



{vnla^r-^Yrs^,{s\{n-s)\r 



Moreover, assuming further the cluster to be in a closed 
packed configuration we can approximate, for A = 0, 
v„ w na"^ I\f2. Notice that for A ^ we would expect 
Vn to change by a constant multiplicative factor which 
would still give the same result for the cluster concentra- 
tions. Clearly a proper analysis of the ri— cluster volume 
would require a MC simulation [7[ . This temperature in- 
dependent approximation gives for n'^ = 100 the results 
shown in Fig. [14] (note that the results have very small 
dependence on n^). 




(w„/cr )" 2"/n! where w„ is the volume of an n— cluster. 



100 



FIG. 14. (Color online) Cluster analysis described in the text 
for ri^ = 100 at various densities. 



Form the figure we can say that our simulation re- 
sults for T* = 0.1 and A = —0.3 have qualitatively 
the same behavior of our oversimplified analysis. This 
justifies the fit of Table U where the Laplace multiplier 
is considered as a fitting parameter. The strong de- 
pendence from the nonadditivity (and on temperature) 
that we observe in the simulation is an indication that 
the approximation of neglecting the excess internal en- 
ergy of a cluster is too severe. One should consider 
z™*™ = e-"^"''/^*(i;„/a3)"-i2"/n!, where f^{T*) = 
L u™(l/a;) dx is the excess free energy per particle of 
the n-cluster and u^{T*) = {(^(t /q^){YJiyj=i (j)ij{rij)) /n 
is the reduced excess internal energy per particle of the 
n— cluster. Note once again that choosing an /™ inde- 
pendent of n would lead to the same oversimplified result 
wc described for the cluster concentrations. What re- 
ally matters is the combined dependence of fn^{T*) on n 
and T* which can be assessed within the MC simulation 
0,|4a|- For example the curves of Figs. [2]and[3]with per- 
colating clusters are better fitted by the following three 
parameters expression xf^ w ;\"+<^" n^'^/nl. 

One thing that can be done is to distinguish amongst 
the clusters of n particles between the ones formed by s 
negative particles and t positive particles with t + s = n, 
as done in Ref. |25] . in order to be able to approximate 
analytically the intra-cluster excess free energy per par- 
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tide 



intra ^ ^ intra 

^n / J s,n—s 



1 



1 



's,t 



xe 



!t! fj3{s+t-i) 



dro . . .dv 



s+t 



(4.6) 



(4.7) 



where the configurational integral goes only over the rel- 
ative positions and it covers the region Vig^t of cluster 
configuration space. This way one can quantitatively [25[ 
estimate how the Tani and Henderson theory [59| devi- 
ates from the exact MC results. 

We immediately see how zl^f'^'^ oc K/a^ becomes big- 
ger and bigger as A — )■ —1 and the same holds for 
all the zljT'l'^"' which clearly dominate over all the oth- 
ers zl^*^"' with s ^ t. This qualitatively explains the 
Fig. 13] as is shown in Fig. [15] where we show the re- 
sults from the approximation described in Appendix |X] 
for n" = 30, T* = 0.1, p* = 0.45, a = 1.5, b = 0.9 and var- 
ious values of A (note that the results have very small 
dependence on n^). 




Ij I I I I I— L. 



100 



FIG. 15. (Color online) Cluster analysis described in the text 
for n" = 30, T* = 0.1, p* = 0.45, a = 1.5,6 = 0.9 at various 
values of A. 



V. CONCLUSIONS 

In conclusion we have performed NVT MC simula- 
tions of the RPM with non-additive hard-sphere diam- 
eters outside the coexistence region with particular em- 
phasis on the clustering properties. In order to establish 
whether the cluster analysis falls outside the gas-liquid 
coexistence region for a given value of the nonadditiv- 
ity we accurately determined the binodals of the non- 
additive fluid using the Gibbs ensemble method after a 
density distribution function analysis to get a first in- 
sights on the shifts of the coexistence region with the 



nonadditivity. It turned out that a negative nonadditiv- 
ity tends to shift to higher temperatures and higher den- 
sities the critical point whether a positive one shifts it to 
lower temperatures and densities. The law of correspond- 
ing states does not seem to be strictly fulfilled over an 
extended region below the critical point for A = 0, ±0.1. 
Our results can be used as a theoretical support to the 
analysis of experimental work on room temperature ionic 
liquids |4l| - |43| where shifts in the liquid-liquid binodals 
akin to ours are observed as a function of the kind of 
solvent used in the ionic mixture. 

From the cluster analysis, we were able to distinguish 
between two kind of behaviors for the cluster concen- 
trations. When we do not observe percolating clusters 
during the simulation, the curves for the cluster concen- 
trations as a function of the cluster size are independent 
of the number of particles used in the simulation. When 
we observe percolation during the simulation the curves 
depend on the number of particles used in the simulation 
but obey a straightforward scaling with A^ relationship. 

At low densities the negative non-additive fluid has 
stronger clustering than the pure RPM whereas at high 
densities is the positive non-additive fluid to have the 
strongest one. The positive non-additive fluid is the first 
one reaching the percolating clusters upon an increase of 
density. This certainly depends on the fact that for a 
positive nonadditivity the ions have less space where to 
move at a given density and, due to he presence of two 
opposite charged species, there is a competition between 
the tendency to clustering driven by the Coulomb inter- 
action and tendency to demixing due to entropic reasons. 
A negative nonadditivity tends to favor the formation of 
the neutrally charged clusters starting with dipolar ones. 
Traces of these features can also be read from an analysis 
of the partial radial distribution function and structure 
factors. Our clustering results can be summarized ob- 
serving that at high density for a positive deviations from 
additivity we have more clustering than in the additive 
model whereas for a negative deviation from additivity 
we have less clustering than in the additive model. At 
low density the reverse behavior is found. These results 
can be explained by the following arguments: at high 
density a positive nonadditivity leaves less effective vol- 
ume to the particles and a negative nonadditivity leaves 
more effective volume relative to the additive model; at 
low density a negative nonadditivity favors the formation 
of neutral clusters and a positive nonadditivity favors the 
competition between the tendency to demixing in a neu- 
tral mixture and the tendency to microscopic intermixing 
of the two species favored by the Coulombic interactions. 
These observations are in agreement with the fact that 
the energy of a cation-anion pair at contact increases for 
positive nonadditivity and decreases for negative nonad- 
ditivity. 

A simple, temperature independent, clustering theory 
where we regard the clusters as forming an ideal gas and 
we approximate the n— cluster as an ideal ensemble of n 
particles in a closed packed configuration can be used to 



15 



qualitatively explain the cluster concentrations observed 
at not to high density and absolute value of the nonaddi- 
tivity. In order to qualitatively explain the prevalence of 
the neutral clusters in the negatively non-additive fluid 
it is necessary to refine the approximation at the intra- 
cluster level. 

In the future it would be desirable to make quantitative 
the comparison between clustering theory and MC exact 
results. It would be also interesting the determination 
of the percolation threshold as a function of nonadditiv- 
ity. In the temperature density phase diagram, one can 
determine the percolation threshold calculating the frac- 
tion of configurations with percolating clusters within the 
NVT simulation as a function of density for two systems 
of different size N . A point of the percolation threshold 
results then from where the curves of the two systems 
meet. 



Appendix A: Approximated intracluster 

configurational partition function for negative 

nonadditivity 



Let us call the anions i_ = 1_, . . . , s_ and the cations 
j+ = 1+, . . . , i-i_. From Eq. (|4.7p in the main text follows 

c/™ = ^^^ / ri ^ri..- n *^>- 

t 

n 

11/" ,', ^ 






-2AbM^,^ -p|- ^+\B/r,^i_ 
t t 



1=2 



k=l 






, + AB/n 3^ 



where wc approximated e ^'^/'^ Ki 1 which is justified at 
high T* < 1/2(1 + A) or low \b- Now we observe that 



for example 'ri^,2_ ~ |ri+i_ +ri_2_ | with ri_2_ > f and 
g+ b/i'i+2_ ~ 2. So that for negative nonadditivity we 
can further approximate 

•"'^t.t 1=2 k=l 



1 f ,*, * 

"'S2t,t ;^2 fc=l 



1 1 

i!2 



x[]e+^- /'••+' 



«(2i)H2*-i)^^^^^^ 



,t!2 



-{K/KoT , 



(A2) 



where the factor (2i)''(^* ^^ takes into account the volume 
of Clt,t with b a free parameter and 

/•Ab/2 /•Ab/2 

K/Ko^ r^e+^^'^'drl r^dr, (A3) 

./acr(l+A) ./acr(l+A) 

with a a second free parameter. With the same approxi- 
mations we can say 



l„ I 4.\b(8+t-l) 

^ntra OC LJLJ_ (T^/A'o)™^"'*^ . (A4) 



's,t 



s!i! 
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